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1  Introduction 


Understanding  and  predicting  the  growth  of  unstable  disturbances  in  the  Earth’s  atmosphere 
is  recognized  by  meteorologists  as  perhaps  the  most  important  factor  in  successful  weather 
prediction.  The  pioneering  work  on  baroclinic  instability  by  Charney(1947)  and  Eady(1949) 
demonstrated  that  the  eigenfunctions  of  simple  linearized  forms  of  the  equations  of  motion 
yield  eigenvalues  which  correspond  surprisingly  well  with  observed  time  and  space  scales 
of  mid-latitude  disturbances.  In  subsequent  years  linear  stability  analysis  of  various  forms 
of  the  linearized  equations  became  the  foundation  for  rapid  advances  in  the  understanding 
of  atmospheric  instability.  Beginning  in  the  1950’s  numerical  computation  made  a  critical 
contribution  to  this  progress  by  making  it  possible  to  solve  eigenvalue/eigenvector  problems 
that  were  intractable  with  analytic  solution  methods. 

All  of  the  early  studies  of  atmospheric  instability  were  based  on  finding  the  normal 
mode  eigenfunctions  of  some  form  of  a  linear  evolution  operator  A  ,  formulated  for  a  steady 
state  mean  flow.  However,  Lorenz  (1965)  showed  that  the  eigenvalues  of  A'^A,  called  the 
singular  value  eigenvalues  of  A,  can  be  much  larger  than  the  normal  mode  eigenvalues  of  A 
itself.  A"^  is  the  transpose,  or  adjoint,  of  A.  The  eigenvectors  associated  with  the  singular 
values  are  the  singular  vectors  of  A.  In  a  meteorological  sense  the  largest  singular  values 
and  their  associated  singular  vectors  are  optimal  in  that  they  maximize  linear  disturbance 
growth  over  a  given  evolution  (prediction)  period.  Also,  because  the  product  operator  A'^A 
is  symmetric,  solutions  for  non-steady  basic  states  are  potentially  more  tractable,  at  least 
in  a  linear  algebra  sense.  However,  the  completely  general  singular  vector  problem,  for  an 
A  with  a  realistic  time  and  space  varying  basic  state  trajectory,  produces  linear  systems 
of  enormous  size  that  require  considerable  computational  effort  and  special  methods  for 
solution.  For  example,  the  evolution  operator  A  for  a  numerical  weather  prediction  model 
of  relatively  modest  resolution  can  easily  yield  linear  systems  of  0(10'^).  Only  recently. 
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therefore,  has  computational  power  made  practical  meteorological  applications  of  singular 
vectors  feasible.  Lacarra  and  Talagrand  (1988)  made  singular  vector  calculations  for  a 
barotropic  shallow  water  model,  and  Molteni  and  Palmer  (1993)  did  similar  calculations 
with  a  barotropic  model  and  a  3-level  quasi-geostrophic  model.  Buizza  et  a/.  (1993)  were  the 
first  to  compute  singular  vectors  for  a  low  resolution  (T21)  numerical  weather  prediction 
(NWP)  model  based  on  the  primitive  equations.  Ehrendorfer  and  Errico  (1994)  made 
similar  calculations  for  a  limited  area  NWP  model. 

Singular  vectors  have  potential  important  application  to  a  wide  variety  of  practical 
meteorological  problems.  In  recent  years  operational  NWP  centers  have  begun  running 
ensemble  prediction  systems,  in  which  each  member  of  the  ensemble  is  an  NWP  model 
integration  run  from  initial  conditions  which  differ  from  those  of  other  ensemble  members 
by  some  suitably  chosen  perturbation  vector.  Ideally,  the  envelope  of  forecast  trajectories 
from  the  ensemble  forecasts  will  occupy  the  full  uncertainty  spectrum  for  this  set  of  numer¬ 
ical  forecasts.  Singular  vectors  are  a  particularly  attractive  source  of  these  perturbations 
because  they  are  orthogonal  to  one  another  and  represent  the  potentially  fastest  growing 
modes  of  instability  for  their  associated  basic  states.  These  instabilities  can  be  expected 
to  give  maximum  divergence  of  forecasts  running  from  initial  conditions  which  differ  by 
combinations  of  singular  vectors. 

Singular  vectors  are  also  important  tools  for  analysis  of  numerical  model  prediction 
error.  It  is  commonly  understood  that  uncertainty  in  the  initial  conditions  (analysis  error) 
contribute  to  the  forecast  errors  that  grow  with  time  in  an  NWP  model.  However,  knowledge 
of  the  geographical  distribution  of  analysis  error  is  not  sufficient  for  understanding  of  why 
these  forecast  errors  grow.  It  is  also  necessary  to  know  where  the  model  is  sensitive  to  errors 
in  its  initial  conditions.  The  geographical  co-location  of  analysis  error  and  forecast  model 
error  sensitivity  is  the  critical  combination  that  leads  to  rapid  growth  of  forecast  error. 
Singular  vectors  associated  with  the  largest  singular  values  are  ideal  representations  of  the 
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forecast  model  error  sensitivity  patterns,  since  each  singular  vector  shows  the  geographical 
distribution  of  perturbation  growth  for  the  associated  singular  value. 

This  report  describes  the  formulation  of  the  singular  vector  system  based  on  the  Navy 
Operational  Global  Atmospheric  Prediction  System  (NOGAPS)  spectral  forecast  model 
(Hogan  and  Rosmond  1991).  In  section  2  a  brief  overview  of  the  NOGAPS  spectral  model 
is  given.  Section  3  describes  the  linearized  version  of  the  NOGAPS  spectral  model  about 
a  time  and  space  varying  basic  state,  often  called  the  tangent-linear  model  (TLM).  Section 

4  describes  the  adjoint  model,  which  is  the  implicit  matrix  transpose  of  the  TLM.  Section 

5  describes  the  procedures  used  for  testing  TLM  and  adjoint  models  and  presents  some 
testing  results.  Section  6  describes  the  formulation  of  the  singular  vector  system,  which 
is  a  linear  operator  based  on  a  combination  of  the  TLM  and  adjoint  operators.  Section  7 
presents  the  results  of  example  singular  vector  calculations.  Section  8  is  a  summary  and 
gives  potential  future  directions  for  singular  vector  and  adjoint  model  applications. 


2  The  NO  GAPS  spectral  model 

The  basis  of  the  singular  vector  system  is  the  non-linear  model  from  which  the  TLM 
and  adjoint  models  are  developed.  We  use  the  NOGAPS  spectral  forecast  model,  described 
in  detail  by  Hogan  et  al.  (1991).  The  model  is  representative  of  the  sophisticated  global 
atmospheric  models  run  by  major  operational  weather  forecast  centers  and  meteorological 
research  groups  around  the  world  for  numerical  weather  prediction,  data  assimilation,  and 
climate  simulation.  NOGAPS  is  the  heart  of  the  Navy’s  operational  weather  prediction 
system  run  at  Fleet  Numerical  Meteorological  and  Oceanographic  Center  (FNMOC).  The 
model  contains  diabatic  parameterizations  of  all  of  the  important  physical  processes  of  the 
atmosphere,  such  as  cumulus  parameterization,  cloud  and  radiation  interactions,  and  the 
planetary  boundary  layer.  Linearizations  of  all  of  these  parameterizations  would  be  desir- 
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able,  but  the  complexity  of  these  schemes  is  daunting,  and  only  a  few  attempts  at  linearizing 
diabatic  processes  for  NWP  models  have  been  made,  e.g.  Vukicevic  and  Errico  (1993).  Fur¬ 
thermore,  mid-latitude  atmospheric  instabilities  are  dominated  by  the  dry  dynamics,  and 
experience  to  date  (Buizza  1994)  has  shown  that  a  linear  model  with  no  diabatic  processes 
except  a  very  simple  planetary  boundary  layer  gives  quite  satisfactory  results. 

For  the  above  reasons  the  linearizations  described  in  the  next  section  will  be  restricted 
to  the  dry  dynamics  with  a  simplified  planetary  boundary  layer  (PBL),  so  in  this  section 
only  the  details  of  the  dynamical  equations  are  presented,  as  formulated  for  the  NOGAPS 
global  spectral  model.  Complete  description  of  the  rest  of  the  model  formulation  is  in 
Hogan  et  a/.  (1991).  Parts  of  the  model  which  are  already  inherently  linear,  such  as  the 
semi-implicit  scheme,  will  not  be  discussed.  Likewise  the  normal  mode  initialization,  whose 
non-linear  operations  are  essentially  identical  to  those  in  the  dry  dynamics,  will  not  be 
discussed.  Note  that  even  though  the  linearization  is  restricted  to  the  dry  dynamics,  the 
basic  state  trajectory  forecast  about  which  TLM  or  adjoint  model  perturbations  vary  is 
always  generated  by  a  full  physics  version  of  the  non-linear  model,  usually  at  relatively  high 
resolution.  This  ensures  that  the  basic  state  trajectory  will  give  the  best  possible  synoptic 
forecasts  and  realistic  linear  instability  growth  behavior. 

2.1  The  basic  equations 

The  NOGAPS  spectral  forecast  model  is  a  primitive  equation  model  formulated  in  spherical 
coordinates  in  the  horizontal  and  a  hybrid  vertical  coordinate  similar  to  that  described  by 
Simmons  and  Striifing  (1981).  The  horizontal  coordinates  are  the  longitude.  A,  and  the 
latitude,  (f.  The  vertical  coordinate  is  normalized  pressure  represented  by  the  variable  r) 
which  ranges  from  0.0  at  the  model  top  to  1.0  at  the  surface.  If  ptop  is  the  pressure  at  the 
top  of  the  model  atmosphere  and  Ps  is  the  terrain  pressure  then  the  pressure  p  is  a  function 
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of  77  given  by 


p  =  ^(77)  +5(77)77, 


(1) 


where  tt  =  ps  —  Ptop-  The  functions  ^(77)  and  B{p)  are  any  two  functions  defined  on  the 
interval  0.0  to  1.0  with  the  boundary  conditions: 


4(0) 

—  Ptop 

4(1) 

~  Ptop 

B(0) 

-  0.0 

B(l) 

=  1.0 

(2) 


where  ptop  is  the  model’s  top  pressure,  ps  the  surface  pressure  and  =  ps  —  ptop- 

The  continuity  equation  for  the  conservation  of  mass  in  this  coordinate  system  is 


(3) 


where  u  is  the  horizontal  velocity  vector.  Integrating  (3)  with  the  top  and  bottom  boundary 
conditions 

77(0)  =  77(1)  =  0,  (4) 


yields  the  tt  tendency  equation 


dn 


(5) 


where  dp  a  function  of  A,  p,  and  77.  We  obtain  the  vertical  motion  equation  by  integrating 
(3)  from  0  to  77  and  substituting  for  the  dTr/dt  term  with  the  right  hand  side  of  (5),  yielding 

("I) -  r  (“I)  ('> 

The  thermodynamic  energy  equation  in  terms  of  potential  temperature  6  is 


dt 


u  dd  V  39 


a  cos  ipd\  a  dp 


.  dp 

% 


dp 


+  Qe-, 


(7) 
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where  a  is  the  Earth’s  radius,  u  and  v  are  the  zonal  and  meridional  wind  components, 
respectively,  and  Qs  is  the  diabatic  forcing  due  to  radiation,  latent  heat  release  processes, 
horizontal  diffusion,  and  vertical  mixing. 


The  moisture  conservation  equation  is 


dt 


u  dq  V  dq 


.  dp 

dp 


V 


(8) 


a  cos  (fdX  a  dip 

where  the  forcing  term  Qq  is  due  to  condensation/evap oration  processes  and  turbulent  and 
cumulus  vertical  mixing. 

We  write  the  hydrostatic  equation  in  the  form: 


^  =  -C0 

dP  '  ’ 

where  4>  is  the  geopotential  and  P  is  the  Exner  function 


(9) 


\Po 


(10) 


In  (10)  po  is  1000  mb  and  k  =  R/cp  is  the  ratio  of  the  gas  constant  R  to  the  heat  capacity 

Cp. 

The  dependent  variables  describing  the  motion  field  in  NOGAPS  are  the  vorticity  C 
and  divergence  D.  This  choice  is  a  routine  one  for  global  spectral  models  because  as  scalar 
quantities  C  and  D  are  easily  expandable  in  terms  of  spherical  harmonics.  To  represent  the 
relationship  between  C,  D  and  the  horizontal  wind  components  it  is  convenient  to  define  the 
operator  a{g,  h),  which  operates  on  any  two  functions  g  and  h,  as 

(11) 


.  1  dg  ^  dh 

a[g,  h)  =  ^ - XTTT  +  7:“! 


\  —  pp  d\  dp 

where  p  =  sin  (p.  The  vorticity  and  the  divergence  are  expressed  as 

(  =  a{V,-U), 


and 


D  =  a(U,V), 


(12) 


(13) 
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where  we  have  defined  the  scaled  wind  velocity  components 


U  =  u 


COSip 


(14) 


and 


V  =  v 


cosp 

a 


Similarly,  the  tendency  equations  for  vorticity  and  divergence  are 


and 


where  the  functions  G,  H,  and  I  are: 


G  =  U{C  +  f)  + 
H  =  V(C  +  /)- 

I  = 


dp 

% 


a 


.  dp 

2  /C/2  +  V2 


'9v\  ^  .  ^  fdp\ 

^dp  J  0?  \d7r  J  \dp  J  a 

'dU\  ,  n 

^dp  j  0?  \dTi)  \d\)  a  ’ 


(15) 


(16) 


(17) 


(18) 

(19) 

(20) 


i-mH  2 

and  /  is  the  coriolis  parameter  and  Qu  and  Qy  are  diabatic  forcing  terms  due  to  surface 
friction  and  vertical  mixing  of  momentum. 


2.2  Vertical  differencing 

NOGAPS  uses  second  order  finite  differencing  for  representing  vertical  derivatives.  The 
distribution  of  variables  on  the  model’s  vertical  grid  is  shown  in  Figure  1.  The  dashed  lines 
denote  ’full’  levels  and  the  solid  lines  ’half’  levels  which  are  at  the  interfaces  between  model 
layers.  The  discrete  analogue  to  (1)  is 

Pk+l/2  =  Afc+i/2  +  Bk+il2TT,  (21) 
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Figure  1:  Terrain  following  vertical  coordinate  for  18  level  NOGAPS  model.  Solid  lines  are 


’half’  levels  , dashed  lines  are  ’full’  levels. 


subject  to  the  boundary  conditions  of  (2).  Model  layer  thicknesses  are  given  by: 


Apk  —  AAk  +  ABkTT 


(22) 


The  discrete  analogue  to  (5)  is 


Cj'K 

=  >  (23) 

1=1 

where  L  is  the  number  of  model  layers.  Analogous  to  (6)  in  the  continuous  case,  the  vertical 
velocity  in  the  hybrid  coordinate  system  is 

■fi^]  =  Bmp  y  (A.4,v  ■  u,  +  AB,V-  (jru,)) 

L  mMP  tl 

k 

-  ^(AAiV-ui  +  ABiV-(7rui))  (24) 

/=i 

At  the  half  pressure  level,  the  Exner  function  P  is: 


p  _  (  Pk-l-l/2 

P^k+l/2  —  I  — - 


and  at  full  levels  is  given  by  Phillips  (1974)  as 


Pk  = 


(k  +  1)p5 


Pt'tl/2-Pt-h2 


From  (21),  (25),  and  (26)  the  term  {dP/dTr)k  in  (18)  and  (19)  becomes,  after  a  little  algebra, 

dPk  _  Pk+i/2  (Pk+i/2  -Pk)  +  Bk-i/2  [Pk  -  Pk-1/2) 

d-K  Apk 

Vertical  advection  requires  that  we  define  half  level  values  of  U,  V,  0  and  q.  These 
choices  are  important  for  energy  conservation.  The  general  form  for  vertical  advection  in 
the  discrete  equations  is 

/ Xk+1/2  -  Afc \ 

V[^5?7j  dp)^  L^^Jjt+i/2V  ^Pk  ) 


+  77 


.  dp 

1 

f  Xk+1/2  - 

dp 

1 

k+1/2 

i,  Apk 

.  dp 

1 

(Xk-Xk- 

^  dp 

1 

fc-1/2 

\  Apk 

where  X  represents  any  dynamic  variable. 

For  the  thermodynamic  energy  equation  the  correct  choice  for  half  level  6  is 


^fc+1/2  =  ^fc+l  +  lVfc+1/2  {Ok  —  ^k+l) 
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where 


(30) 


.A,  ^Pk+l/2-Pk 

m+1/2  -  -B - 

(Haltiner  and  Williams  1980).  This  ensures  total  energy  conservation  in  conversions  between 
kinetic  and  potential  energy  in  the  model  atmosphere.  Note  that  this  is  not  a  linear  inter¬ 
polation  in  P.  Details  can  be  found  in  Hogan  et  a/.  (1991).  Using  (29)  the  discrete  form  of 
(7)  is 


dt 


=  dkD,-a{Ukek,Vkek) 


.  dp 

% 


'6 


'k+l/2 


-e. 


Jfc+l/2 


.  dp 


'6h  —  Oh- 


'fc-l/2 


ik-lj2 


^Pk 


+ 


(31) 


where  we  have  used  the  identity  u-V^  =  V-(u6')  -  0V-u  and  a  generalization  of  (13)  to 
redefine  the  horizontal  temperature  advection  terms. 

For  half  level  U  and  V  the  linear  interpolation 


-^k+l/2 


(32) 


where  X  is  one  of  the  two  scaled  wind  components,  guarantees  kinetic  energy  conservation 
for  vertical  advection  of  momentum.  This  is  equivalent  to  W  =  0.5  in  (30)  .  The  same 
form  is  used  to  define  half  level  g,  although  this  can  lead  to  problems  because  it  cannot 
guarantee  that  positive  definite  moisture  is  preserved  by  the  vertical  advection.  However, 
because  we  are  mainly  concerned  with  the  linearization  of  the  dry  dynamics  of  NOGAPS, 
we  will  not  concern  ourselves  with  this  problem  here.  Using  (32)  for  defining  half  level  U, 
V,  and  q,  the  discrete  forms  of  (16),  (17),  (18),  (19),  (20),  and  (8)  are,  respectively, 

=  -a{Gk,Hk),  (33) 

=  a{Hk,-Gk)-V^<Pk  +  h)-  (34) 


dCk 

dt 

dPk 

dt 


Gk 


Ufc(a  +  /)  + 


‘•^1  (Vk^i-Vk\ 

.'^^^Jfc+i/2V  2Apfc  ) 


k-l/2 


\  2Apfc  J 
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<f>k-^s  =  Cp6k{Pk+i/2  -  Pk),  k  =  L 

4>k  ~  (f’k+l  =  Cp^k+l/2{Pk+l  ~  -Pfc))  L  —  1  >  k  >  1 
where  (ps  is  the  terrain  geopotential  and  6k+i/2  is  the  half  level  potential  temperature  given 
by  (29). 


2.3  Spectral  transforms 

The  NOGAPS  model  equations  presented  above  are  solved  in  a  sequence  of  steps  which 
require  that  model  variables  be  represented  both  in  ’spectral  space’  and  ’grid  point  space’. 
Spherical  harmonic  transforms  are  the  natural  choice  for  the  equations  expressed  in  spher¬ 
ical  coordinates.  Spherical  harmonics  are  the  associated  Legendre  polynomials,  P^  (p) 
multiplied  by  the  complex  Fourier  series,  The  subscript  n  is  the  total  wavenumber 

and  superscript  m  the  zonal  wavenumber.  The  series  are  truncated  assuming  a  triangular 
truncation  with  M  indicating  the  total  number  of  resolvable  waves.  If  is  a  set  of 

complex  spherical  harmonic  coefficients  for  a  variable  X,  then  the  spectral  to  grid  spectral 
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expansion  of  X  is 


M  \  M  1 

X{\,,n<ril.,t)=  E  E  5.”  <=■”"'.  (40) 

m=~M  ji^\m\ 

where  X  is  any  of  the  model  dependent  variables  tt,  D,  6,  or  q  at  any  time  t  The 
subscripts  I  and  j  represent  the  discrete  grid  longitudes  and  latitudes,  respectively,  of  a 
suitably  chosen  transform  grid,  and  the  rjk  are  discrete  vertical  coordinate  values.  The  inner 
sum  in  (40)  yields  a  set  of  complex  Fourier  coefficients,  so  the  outer  summation  is  a  discrete 
backward  Fourier  transform  which  is  solved  by  fast  Fourier  transform  (FFT).  The  choice  of 
the  transform  grid  resolution  associated  with  a  particular  value  of  M  is  an  important  one 
to  ensure  accuracy  and  integral  property  conservation.  Details  of  this  choice  are  in  Hogan 


et  a/.(1991). 

The  grid  point  to  spectral  transform  is  obtained  from  (40)  by  applying  the  orthogon¬ 
ality  property  of  the  spherical  harmonics: 

Sr  (l.  <)  =  ^  /-,  [C  ^  (41) 

The  expression  inside  the  square  brackets  is  a  forward  Fourier  transform  also  solved  by 
FFT.  After  the  FFT,  (41)  is 

('). «)  =  ^  1^  (i*.  1)1  ('i)  (^2) 

Where  [X  (/x,  r],  t)]  are  the  coefficients  from  the  Fourier  transform  described  above.  The 
integration  in  (42)  is  performed  using  the  method  of  Gaussian  quadrature,  for  which,  given 
any  polynomial  ^'(/x),  of  degree  2J  —  1  or  less,  the  integral  of  g  from  -1  to  1  is  computed 
exactly  as 

/•i  ^ 

“'-i  j=\ 

where  iij  are  the  so-called  Gaussian  latitudes  and  Wj  the  Gaussian  quadrature  weights. 
With  this  the  discrete  form  of  (42)  is 

(3M-M)/2 

j=i 
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The  upper  limit  on  the  summation  is  the  number  of  Gaussian  latitudes  necessary  to  ensure 
that  the  integral  in  (42)  is  evaluated  exactly. 

Horizontal  derivatives  of  spectral  variables  are  found  by  differentiating  (40),  which 


yields 


and 


m=~M 


M 


7n=—M  [_n=|m| 


dn 


r,imXi 


^imXi 


(44) 


(46) 


It  is  also  necessary  to  have  transform  expressions  relating  the  wind  field  to  vorticity 
and  divergence.  The  scaled  wind  components  U  and  V  are  defined  in  terms  of  the  stream 
function  ip  and  velocity  potential  x  ^ 


I  dx  \  —  y?  dip 

dX  oX  dpi  ’ 


(46) 


and 


1  -  dip 


(47) 


a?  dp,  a?  dX 

and  D  are  related  to  ip  and  X  by  C  =  and  D  =  V^x-  Spherical  harmonics  are 
eigenfunctions  of  the  Laplacian  operator  V^,  with  eigenvalues  — n(n  +  l)/a^.  Therefore 

n(n  +  1) 


c 


m 

n 


Dl 


n{n  +  1) 


i/;”" 
'rn  ? 


m 
An  ’ 


Combining  (46)  and  (47)  with  (48)  and  (49)  gives  the  spectral  expansions 


M 


U{Xi,pj,r]k,t)  = 

m=~M 


M  1 

y  — = — 

I  I  I  nin  +  1) 


Br(>?*,i)p„”(«) 


^im\i 


M 


+  (i-M^)  E 

m=~M 


M 

E 


1 


^imXi 


(48) 

(49) 


(50) 
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and 


'3 )  ^imXi 


m=—M  ji=\m\  ^  ^ 

M  [  M  1  WP^/'//'ll 

-  (1-.^)  E  (5» 

m—~M  [_n=|m|  v  '  /  j 

The  spectral  coefficients  of  the  vorticity  and  the  divergence  can  be  obtained  from  the  grid 
point  fields  of  U  and  V  by  using  the  spectral  representations,  which  are  given  by  (50)  and 
(51),  the  orthogonality  of  the  expansion  functions,  and  the  zero  boundary  conditions  of  U 
and  V  at  the  poles.  Following  the  procedure  outlined  in  (42)  and  (43)  this  yields 

{3M+l)/2  z'  P^fn) 

C(Vk,t)  =  E  Wj\imT^[V{iij,r]k,t)]  "  V 


+  [U{iJ,j,r]k,t)] 


and 

(3M+l)/2  r  P^(u) 

;  =  1  I  ^  ^3 

-  (53) 

where  as  before  the  T  represents  the  complex  Fourier  coefficients  of  the  respective  wind 
components  in  the  square  brackets. 


2.4  Time  integration 

The  NOGAPS  spectral  model  uses  leapfrog  time  differencing  with  a  Robert  time  filter.  This 
is  particularly  important  for  the  adjoint  model,  which  must  integrate  backwards  in  time  in 
a  way  that  is  consistent  with  the  forward  time  integration  scheme.  An  integration  is  started 
with  a  forward  time  step: 

r)Y 

=  X,^  +  dt  ^  ,  (54) 

.  ^  J  t=0 
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followed  by  leapfrog  and  time  filter  steps, 


-^t=2 


Xtz=o  d"  ^dt 


dX 


dt 


j  t=i 


Xt-i  +  e  [Xt=2  ~  2X(=i  +  Xt=Q 


(55) 

(56) 


then 


= 

= 

and  so  on.  dt  is  the  time  step  and  e  is  the  time  filter  coefficient:  0  <  e  <  0.25.  The 
are  the  time  filtered  values  that  replace  the  existing  X  (t) .  Distinguishing  between  the 
and  X(t)  is  an  important  detail  for  the  adjoint  of  the  time  integration  scheme. 


Xt=i  +  2dt  ^ 

[9tu=2 

Xt=2  +  ^  [Xt^s  —  2Xt=2  +  Xt-i] 


(57) 

(58) 

X(i) 

X(t) 


3  The  tangent  linear  model 

The  linearization  of  the  set  of  equations  given  in  section  2  above  is  straightforward.  The  only 
novel  feature  of  the  resulting  linear  system  of  equations  is  that  they  describe  the  behavior  of 
perturbations  about  a  basic  state  which  varies  in  time  and  all  three  space  dimensions.  This 
means  that  the  coefficients  of  the  linear  model  that  depend  on  this  basic  state  will  be  tangent 
to  the  nonlinear  model  trajectory  in  phase  space,  hence  the  term  tangent  linear  model,  or 
TLM.  The  TLM  is  nearly  always  an  extremely  close  approximation  to  the  nonlinear  model 
unless  the  linear  approximation  breaks  down  even  for  very  short  time  and  space  scales. 

A  very  important  practical  implication  of  the  time  varying  basic  state  is  that  during  a 
trajectory  generating  non-linear  model  integration  the  forecast  variables  (coefficients)  must 
be  writing  to  a  file  at  sufficiently  high  frequency  to  resolve  the  trajectory.  This  is  often  every 
time  step  of  the  non-linear  model  integration.  Then  during  a  TLM  model  (or  adjoint  model) 
integration  based  on  this  trajectory  the  file  must  be  read  at  this  frequency  to  retrieve  the 
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coefficients  in  the  appropriate  order.  This  implies  a  tremendous  I/O  bandwidth  overhead 
for  TLM  and  adjoint  model  integrations  far  in  excess  of  normal  NWP  model  integrations. 
During  singular  vector  calculations  this  I/O  requirement  can  be  extremely  burdensome  and 
often  requires  alternate  strategies  such  as  storing  the  trajectory  in  memory,  if  possible,  or 
using  special  I/O  devices  with  much  better  performance  than  typical  rotating  mass  storage 
devices.  This  extra  burden  on  TLM  and  adjoint  integrations  led  Errico  et  al.  (1993)  to 
discuss  the  impact  of  less  frequent  trajectory  update  on  TLM  integrations. 

To  demonstrate  a  simple  linearization  example,  consider  the  simple  quadratically  non¬ 
linear  equation  y  —  ab.  Assume  a,  b,  and  y  are  composed  of  basic  state  components  a,  6,  y 
and  perturbations  a' ,  b',  and  y'  about  the  basic  state.  We  have 

y  +  y'  =  {a  +  a')(b  +  b'). 

Collecting  terms  and  assuming  y  =  ab  yields 

y'  =  ab'  -f  bo!  +  a'b' . 

The  linearization  is  completed  by  neglecting  the  a'b'  term,  which  should  be  small  compared 
to  the  linear  terms  if  perturbations  are  small.  For  higher  order  nonlinearities  there  would 
be  additional  terms  to  be  neglected,  but  the  assumption  of  small  perturbations  relative  to 
the  basic  state  variables  is  still  the  basis  of  successful  linearization. 

We  shall  derive  tangent  linear  expansions  of  each  of  the  important  governing  equations 
presented  in  section  2(b)  above.  Intermediate  variables  will  be  defined  where  convenient 
and  to  assist  clarity.  An  important  issue  in  TLM  and  adjoint  model  development  is  coding 
order,  since  the  adjoint  model  is  a  series  of  matrix  transposes  which  are  equivalent  to  exact 
reversal  of  the  logical  flow  of  expressions  in  the  TLM  model.  We  will  therefore  present 
equations  in  essentially  the  same  order  they  occur  in  the  TLM  model  code.  The  complete 
TLM  code  is  available  upon  request  from  the  author.  We  will  follow  the  normal  convention 
of  representing  a  perturbation  variable  as  (  )^  and  a  basic  state  variable  as  (  )■  A  nonlinear 
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NOGAPS  model  variable  expansion  in  discrete  form  is  therefore 

where  l,j,k  are  the  indices  for  longitude,  latitude,  and  vertical  coordinate,  respectively. 
Note  the  full  space  and  time  dependence  of  the  basic  state.  For  convenience  in  the  following 
we  will  drop  all  but  the  vertical  index. 


3.1  Spectral  to  real  transforms 

The  time  stepping  process  in  the  NOGAPS  spectral  model  is  done  with  the  variables  defined 
in  ’spectral  space’.  The  same  is  true  for  the  NOGAPS  TLM  and  adjoint  models.  Therefore 
at  the  beginning  of  each  TLM  time  step  spectral  variables  must  be  transformed  from  spectral 
to  ’grid  point  space’  to  allow  calculation  of  linearized  grid  point  tendencies.  To  get  grid 
point  TT,  C,  D,  9,  and  q  we  use  (40).  Horizontal  gradients  of  dependent  variables  are  found 
using  (44)  and  (45).  For  the  vector  wind  components  U  and  V  (50)  and  (51)  are  used.  Since 
these  are  linear  transforms  the  TLM  model  needs  no  changes  from  the  nonlinear  NOGAPS 
code.  The  basic  state  variables  are  also  carried  in  spectral  form  so  equivalent  transforms 
must  be  made  to  get  grid  point  basic  state  quantities.  With  these  three-dimensional  grid 
point  fields  for  the  dependent  variables  we  can  compute  all  necessary  diagnostic  quantities 
and  TLM  tendencies. 


3.2  Terrain  pressure  dependent  linearizations 

The  perturbation  forms  of  (21)  and  (22)  are 

Pk+l/2  ~  ^k+112'^  ) 
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and 


Apfc  =  ABfcTr'. 


(60) 


Note  that  for  linear  equations  such  as  these  the  basic  state  does  not  appear.  The  surface 
pressure  tendency  (23)  yields 


^TT 


U  H 


k=\ 


,  (  Uk  dll'  —  dir'  , 


+ 


1  dw^^,  dTf,„ 


(61) 


and  (24) 


.  dp 

% 


1  / 


k  r 


=  -E 

fc+l/2  k=l 


AAk-D'k  +  AB); 


U  k  dir' 


9^'  ,  n 


1  —  aA  op 


1  dW  ,  dW  ,  _ 

+  - 2 

I  -  p^  oX  op  ^ 


dll' 

—  'Dfc+1/2"^- 


dt 


(62) 


Note  the  general  characteristic  of  the  linearized  equations  to  be  more  complex  than  their 
nonlinear  counterparts,  typically  with  each  nonlinear  term  becoming  two  or  more  linear 
terms.  Therefore  the  TLM  is  significantly  more  expensive  to  run  that  the  dry  version  of  the 
nonlinear  NOGAPS. 

Linearization  of  the  transcendental  Exner  functions  (25)  and  (26)  is  most  convenient 
with  a  first  order  Taylor  series  expansion.  Because  P  is  only  a  function  of  the  terrain 
pressure  tt,  we  have  _ 

F  = 


diT  J 


so  the  half  level  perturbation  Exner  function  is 


I^Bk+i/2Pk+l/2 


(63) 


At  full  levels  using  (27)  yields 


P' 

^k 


' Bk+ll2  (Pfc+l/2  —  -Pfc)  +  Bk-i/2  (-Pfc  Pk-1/2) 

^Pk 


TT  . 


(64) 
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To  linearize  (27)  we  use  (60),  (63),  and  (64): 


m 

dn 


Bh+l/2  (■Pfc+l/2  “  +  -Sfc-l/2  [Pk  Pk-l/2j 

^Pk 

Bk+l/2  [P fc+1/2  —  Pk)  +  Pk-1/2  (^Pk  —  P fc-l/2) 


^pI 


ASfcTr'. 


(65) 


3.3  The  hydrostatic  equation 

The  linearization  of  (30)  is 

WUl/2  =  W,+i/2 

which  combined  with  (29)  leads  to 


P' 

^k+ll2 


PL 


p'  —  p' 


L^/c+1/2  -  Pk  Pk^l  -  Pk\ 


(66) 


^fc+1/2  ~  ^'k+l  +  ^fc+1/2  ^  ^fc+l)  +  Wfc+1/2  ~~  )  •  (67) 

The  hydrostatic  equation  (39)  is  linearized  using  the  definition  for  ^^+1/2,  Pk+1/2^ 

P^.  given  in  (67),  (63),  and  (64).  We  have 

^'k  ~  Cp9i^{Pk+l/2  ~  Pk)  +  Cp9k{Pf.^i/2  ~  Pk}i  k  —  L  ^ 

4>'k  ~  ^'k+1  ~  ^P^k+l/2iPk+l  ~  Pk)  +  Cp9k+i/2{Pk-\-l  ~  Pk)^  L  —  1  >  A:  >  1 
The  discrete  hydrostatic  equation  is  an  integral  of  atmospheric  mass  from  the  Earth’s  surface 
upward.  (68)  shows  the  dependence  of  a  layer  k  on  layers  below  it,  e.g.  A:  4- 1.  The  ordering 
of  the  k  index  is  also  explicitly  shown.  Reversing  this  order  will  be  critical  for  the  adjoint 
of  the  TLM  hydrostatic  equation. 


3.4  The  potential  temperature  and  moisture  equations 

To  linearize  any  of  the  multi-level  prognostic  equations  such  as  (31)  we  need  a  linearization 
of  the  vertical  advection  operator  (28).  Using  (60)  and  (62)  the  perturbation  form  of  (28)  is 


.dp]  dX\' 

\^dr]\  dp) I 


.  dp 


J  fc+1/2 


I  -^^+1/2  ~  A'fe 

V  ^Pk 


+ 


.  dp 

% 


J  fc-l/2' 


'Xk-Xk-i/T 
.  ^Pk 


19 


,  dp 
^  dp 


fc+l/2 


+ 


.  dp 
dp 


V 


f  Xk+l/2  —  Xk 
\  {^Pk?  , 

(xUm2-K^ 


TT  — 


.  dp 

% 


'X,-X 


k-l/2 


J  A-1/2' 


(Xpk)^ 


iTT 


J  fc+1/2 


I  ^ 


+ 


dp 

% 


k-l/2^ 


^  y'  __ 

^k-l/2 

.  ^Pk 


(69) 


With  (69)  the  perturbation  forms  of  (31)  and  (38)  are 

^  =  DA  +  »tD',  -  a  {U'A,  V[ik)  -  a  (uAt.VA)  - 


and 


^  =  DA  +  -  a  {U'A,  VA)  -  a  (Dki,  VA,) 


.  dp 

^¥p 


.  dp 

% 


dp, 


+  <3i..  (70) 


1^)  +0;.. 

k 


dp 


(71) 


Because  we  only  has  a  dry  TLM  at  this  point,  the  Qg^  and  terms  only  contain  a  simple 
vertical  mixing  process  in  the  planetary  boundary  layer  (PBL).  No  perturbation  surface 
sensible  or  latent  heat  fluxes  are  included.  Future  applications  will  require  the  inclusion 
of  these  fluxes  as  well  as  some  linearized  moist  physics  parameterization.  The  simplifled 
linear  mixing  parameterization  is  described  in  section  (f)  below. 


3.5  The  vorticity  and  divergence  equations 

The  perturbation  forms  of  the  C  and  D  tendency  equations  (33)  and  (34)  are  derived  by 
first  linearizing  the  expressions  for  Gk  (35)  (36),  and  Ik  (68).  Using  (69)  and  (65),  we 

have 


G'.  =  U'k{C,k  +  f)^UkCk  + 


•^1 

^dv\  dp)^ 


.  (9P\  dw  ,  „  /-apy  M 

Bnji,dp 


d-K  I  .dp 


+  e 


^dTrJkdp 


Q\ 


Vk 


COS^ 

) 

a 


H'k  =  K(a  +  /)  +  na- 


.  dp 

["^^77]  dp)^ 


(72) 
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^,(dP\dTf  -  fdP\' dn  - 


[dTTj^dX 


,  C0S(P 

+  ^C/fc  ^  . 


4 


1-  fx- 


iU'Uk  +  V'Vk). 


With  these  definitions  the  TLM  tendency  equations  for  C  and  D  are 


dD 


dt 


i  =  Q  (ff;. -Gy  -  (0i  + /;) , 


(73) 

(74) 

(75) 

(76) 


The  diabatic  forcing  terms  Qy^  and  in  G'^  and  H^,  respectively,  could  potentially 
include  perturbation  momentum  fiuxes  by  model  parameterizations  such  as  gravity  wave 
drag  and  cumulus  friction.  Currently,  however,  the  terms  only  include  a  simplified  planetary 
boundary  layer  parameterization  that  will  give  us  a  perturbation  surface  drag  and  vertical 
mixing  of  momentum  contribution  to  prevent  the  growth  of  spurious  unstable  disturbances 
in  the  lowest  model  levels  (Buizza  1994).  The  form  of  these  forcing  terms  is  given  in  section 
(f)  below. 


3.6  Diabatic  forcing  terms 

The  surface  flux  and  vertical  mixing  parameterizations  in  the  NOGAPS  spectral  model  are 
based  on  Louis  (1979).  Details  are  in  Hogan  et  al. (1991).  In  the  linearization  of  this  scheme 
we  will  neglect  buoyancy  effects  and  make  other  simplifying  assumptions  described  below. 
The  vertical  diffusion  of  any  variable  X  =  9,  q,  U,  or  V  is: 

dt  ~  ^  dp' 

where  g  is  the  acceleration  of  gravity  and  P  is  the  vertical  flux  of  X.  We  choose  the  following 
TLM  discrete  form: 

^X'f,  _  J^k+ll2  -  ■^fc-l/2 
dt  ^  Apk  ' 
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where  we  have  neglected  the  effects  of  pressure  perturbations,  i.e.,  tt',  on  the  perturbation 
mixing.  Above  the  surface  the  half  level  perturbation  fluxes  JA^+1/2  are  given  by 

^ fc+1/2  “  4+1/2  Pfc+1/2  '^^*/a:+1/2 


K  - 

Aifc+1/2  J 


1  <  A;  <  L-  1, 


(78) 


where  p  is  air  density,  tZ*  is  the  friction  velocity,  ho  is  the  PBL  depth,  :z  is  the  height  above 
the  surface,  and  Az  is  model  layer  thickness.  Z^i/2  ^he  mixing  length  with  an  asymptotic 
limit  A  given  by 


jm  _  kyZk+ll2 

fc+1/2  l  +  kyZk+1/2/^ 

where  ky  =  0.4  is  the  Von  Karman  constant,  /jt+1/2  is  the  similarity  function 


(79) 


/ic+1/2  —  exp 


-Z^fc+l/2\ 

ho  )’ 


(80) 


Perturbation  surface  fluxes  are  given  by 


_  ^Ps  •a*  y/ 
ln(zL/zo) 


(81) 


where  zq  is  the  surface  roughness  length.  The  surface  momentum  flux  perturbations 

for  X  =  U  and  X  =  V.  We  are  not  allowing  perturbation  sensible  and  latent  heat  fluxes  in 
the  TLM,  so  the  X'g  are  zero  (or  X  —  6  and  X  =  q. 

There  are  several  free  length  scale  parameters  in  (78),  (79),  (80),  and  (81).  In  the 
NOGAPS  TLM  we  use  the  values  proposed  by  (Buizza  1994). 


ho  =  1000.0  m  , 

A  =  80.0  m  , 

zo  =  0.05m  (land),  Zo  =  0.0005m  (ocean). 

For  simplicity  we  use  the  same  value  of  A  for  both  heat  and  momentum  mixing.  A  logical 
choice  for  the  friction  velocity  u*  would  be  a  basic  state  u*  generated  by  the  nonlinear 
NOGAPS  model  and  included  as  part  of  the  time  and  space  dependent  trajectory  which  the 


22 


TLM  and  adjoint  models  need.  Buizza  (1994)  proposed  a  simpler  alternative  that  eliminates 
the  need  for  any  trajectory  dependence.  He  set 

u*  =  0.5  m  sec“^  (land),  n*  =  0.2  m  sec“^  (ocean). 


We  have  experimented  with  both  approaches  and  find  little  difference  in  TLM/adjoint  res¬ 
ults.  For  simplicity  reasons  most  of  our  results  to  date  have  been  with  the  Buizza  values, 
but  the  trajectory  method  is  available  as  an  option  in  the  TLM  and  adjoint  codes. 

To  solve  (77),  we  use  a  backward  implicit  leapfrog  time  differencing  scheme.  This 
leads  to  a  simple  adjustment  scheme  for  the  dependent  variables,  as  described  in  Hogan 
et  ol.(1991).  Substituting  (78)  and  (81)  into  (77),  rearranging  terms,  and  dropping  the 
primes  on  the  X’s  yields 

+  i>kXk  -f  CfcWfc+i  =  Xk,  (82) 

where  Xk  is  the  vertical  profile  of  X  before  vertical  mixing,  and 


ai  = 
dk  = 
Ck 


—2dt  g 
—2dtg 


0, 

ApfcA2j,_i/2 

^r+1/2  Pk+1/2  fk+1/2 

Apt  Azt+1/2 


Cl  = 
bi  = 
bk  = 
bL  = 


0, 


T/2  Pi/2  “*  /1/2' 
ApiAzi/2  I  ’ 


1  —  Cl  -I-  2dt  g 

1  dk  Cki 

l-aL-\-2dtg  [ap^^V/zo)]  ■ 


2<k<L, 

l<k<L-l, 


2<k<L-l, 


(83) 

(84) 

(85) 

(86) 

(87) 

(88) 

(89) 


Equation  (82)  is  a  tri-diagonal  system  of  linear  equations  easily  solved  by  Gaussian 
elimination.  The  set  of  coefficients  (83)  -  (89)  is  appropriate  for  the  vertical  fluxes  of  U 
and  V.  For  9  and  q  vertical  mixing  the  same  set  is  used,  except  that  currently  we  have  no 
perturbation  surface  fluxes  of  heat  or  moisture,  so  &£,  =  1  —  for  these  TLM  variables. 
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3.7  Real  to  spectral  transforms  and  time  stepping 

The  preceding  TLM  development  yields  the  linear  tendencies  of  tt,  C,  D,  9,  and  q  in  ’grid 
point  space’.  These  grid  point  tendencies  are  transformed  to  ’spectral  space’  using  (42), 
(52),  and  (53).  The  spectral  tendencies  are  adjusted  for  the  semi-implicit  corrections  as  in 
Hogan  et  a/. (1991).  Finally  a  time  step  is  done  in  spectral  space  as  in  (54)  to  start  a  model 
integration,  or  as  in  (55)  and  (56)  for  a  leapfrog  and  time  filter  step,  and  the  TLM  model 
variables  are  ready  for  the  next  time  step. 


4  The  adjoint  model 

We  define  the  adjoint  of  the  tangent  linear  model  operator  A  as  .  We  also  define  A  as 
a  sequence  of  linear  operations  with  two  levels  of  abstraction.  At  the  first  level  it  is  the 
product  of  matrices  that  represent  the  time  evolution  nature  of  a  model  integration,  i.e., 

A  =  n  A,  (90) 

*:=! 

where  each  Ak  is  the  matrix  operator  at  time  step  ^  in  a  time  integration  of  n  steps.  At  the 
second  level  each  A^  is  broken  down  into  the  product  of  linear  operators,  i.e., 

m 

Afc  =  0/  (91) 

i=i 

where  the  ai  are  the  individual  fortran  DO  loops  and  lines  of  code  of  the  TLM,  here  shown 
to  number  m.  From  the  formal  definition  of  the  adjoint  as  the  matrix  transpose  of  A  we 
have 

Al=U  af,  (92) 

l=m 

as  the  adjoint  of  (91)  and 

A^  =  ri  Al,  (93) 

k=n 


24 


as  the  adjoint  of  (90).  Notice  the  exact  reversal  of  matrix  product  order  between  the 
TLM  and  adjoint  operations.  In  the  following  sections  we  will  preserve  this  ordering  in  the 
discussion  of  the  adjoint  code  development.  The  special  case  of  the  self-adjoint  spectral 
transforms  will  be  discussed  separately  after  the  rest  of  the  adjoint  model  is  presented. 

There  are  two  approaches  to  developing  adjoint  codes  at  the  practical  level,  by  manual 
methods  and  with  automatic  adjoint  generating  software.  Each  approach  has  its  strengths 
and  weaknesses,  and  a  strategic  combination  of  the  two  was  used  for  the  NOGAPS  adjoint 
code  development,  as  discussed  below. 


4.1  Time  step  and  time  filter  adjoint 


In  the  manual  adjoint  coding  method  each  DO  loop  of  the  TLM  is  written  as  a  matrix  oper¬ 
ation,  the  matrix  is  transposed,  and  an  adjoint  DO  loop  is  written  based  on  the  transpose 
matrix.  As  an  example  consider  the  leapfrog  time  step  and  time  filter  of  (55)  and  (56). 
In  practice,  of  course,  this  is  not  how  it  would  be  programmed  in  a  numerical  model.  In 
NOGAPS  we  have 

X'  =  +  2dtX\  (94) 

=  A"  +  e(A'  -  2A”  +  X°),  (95) 

A"  =  A',  (96) 


where  A'  is  the  tendency,  dX/dt,  overwritten  with  the  ’new’  value  of  A,  X°  is  the  ’old’  time 
level,  and  A”  is  the  ’now’  time  level.  In  matrix  form  (94)  expands  to 


- 1 

o 

_ 1 

1 

O 

_ 1 

1 - 

O 

1 

1 

2dt 

1 - 

1 _ 

(97) 
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(95)  becomes 


X' 

1  0  0 

X' 

0  10 

X” 

e  1  —  2e  e 

X'’ 

and  (96)  becomes 


- 1 

_ 1 

1 

1 

o 

1 

_ J 

X” 

1 

o 

\ _ 

X" 

Following  the  steps  outlined  above  for  adjoint  code  ordering,  we  can  write 


(98) 


(99) 


1 _ 

1 

1 

1 _ 

X2  _ 

0 

0 

- 1 

_ j 

X'a 

1 

0  e 

X" 

— 

0 

1  1  -  2e 

0 

0  e 

- 1 

(100) 


(101) 


1  1 
0  2dt 


XL 


(102) 


where  the  Xa  are  adjoint  variables.  Expanding  back  to  equation  form  yields  the  adjoint  of 


(96) 


the  adjoint  of  (95) 


and  the  adjoint  of  (94) 


II 

X'^^x: 

x:  = 

0 

II 

(1  -  2£)X- 

= 

= 

II 

2dtX'^ 

(103) 


(104) 


(105) 
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Now  for  some  observations.  Notice  that  in  (97),  (98),  and  (99)  each  left-hand  side 
variable  was  assigned  to  the  bottom  row  of  its  matrix  expansion.  This  is  required  to  ensure 
proper  ordering  of  the  adjoint  expressions.  Also  notice  that  in  (103),  (104),  and  (105)  the 
first  equation  is  an  accumulation  into  a  variable,  i.e.,  a  variable  occurs  on  both  sides  of  the 
equal  sign.  Therefore  care  must  be  taken  to  ensure  that  this  variable  is  initialized  properly 
at  the  beginning  of  each  adjoint  time  step.  With  almost  no  exceptions  adjoint  variables 
should  be  set  to  zero  at  the  beginning  of  each  time  step,  and  often  at  intermediate  points 
in  the  code  as  well.  The  general  rule  is  that  when  in  the  TLM  a  variable  is  assigned  a 
new  value  with  no  dependence  on  a  previous  value,  then  at  the  corresponding  point  in  the 
adjoint  code  the  companion  adjoint  variable  must  be  set  to  zero.  For  example,  in  (103)  the 
right-hand  side  X'  must  be  set  to  zero  for  correct  results  from  this  sequence  of  equations. 
Failure  to  set  adjoint  variables  to  zero  at  appropriate  times  is  the  overwhelming  source  of 
errors  in  adjoint  models,  as  this  author  can  wearily  attest. 

Clearly  the  preceding  exercise  demonstrates  the  mechanical  nature  of  adjoint  code 
generation,  which  leads  to  the  obvious  question:  what  about  automatic  adjoint  generation 
software?  In  recent  years  considerable  effort  has  be  devoted  to  such  software  development, 
and  it  has  made  adjoint  model  development  much  less  arduous  by  relieving  modelers  of 
the  tedium  of  hand  coding  the  operations  such  as  described  above.  For  the  NOGAPS 
adjoint  model  the  software  of  Giering  and  Kaminski  (1996)  was  used  on  nearly  all  of  the 
TLM  code  described  in  part  3.  Adjoint  software  is  particularly  well  suited  for  operating 
on  individual  subroutines  at  the  bottom  of  a  calling  tree  with  all  communication  to  the 
rest  of  the  model  through  calling  arguments,  i.e.,  no  COMMON.  In  general  the  software 
cannot  reliably  produce  correct  adjoint  code  for  routines  which  contain  subroutine  calls  or 
COMMON  variables.  The  best  strategy  for  adjoint  model  development  is  therefore  to  use 
the  automatic  software  for  clean  code  areas  of  the  TLM,  and  use  manual  methods  in  top 
level  routines  which  contain  inter-procedural  communication.  The  fortran  code  produced 
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by  the  automatic  software  is  also  cumbersome  and  difficult  to  read,  so  hand  editing  is  often 
necessary  to  get  more  esthetically  pleasing  code.  The  automatic  software  is  also  unable  to 
recognize  self-adjoint  constructs  such  as  the  spectral  transforms  in  the  NOGAPS  TLM,  so 
here  too  manual  methods  must  be  used. 

We  will  now  proceed  step  by  step  though  an  adjoint  model  time  step  in  the  appropriate 
order.  We  will  drop  the  subscript  a  for  convenience,  since  in  the  subsequent  sections  ALL 
references  to  dependent  variables  will  be  to  the  adjoint  variables.  For  sake  of  brevity  not 
all  the  details  will  be  presented.  In  particular  the  adjoints  of  the  semi-implicit  algorithms 
and  horizontal  diffusion  will  not  be  shown,  since  they  are  linear  operations  and  are  easily 
handled  by  automatic  adjoint  software.  Each  adjoint  model  step  begins  with  the  adjoint 
of  the  time  step  and  time  filter,  which  is  already  described  above.  Next  is  the  adjoint 
of  the  grid-to-spectral  transform.  We  shall  postpone  the  discussion  of  this  until  the  end 
of  this  section,  since  it  is  easier  to  understand  the  self-adjoint  properties  when  both  grid- 
to-spectral  and  spectral-to-grid  transform  adjoints  are  described  together. 

4.2  Diabatic  forcing  adjoint 

The  simple  form  for  TLM  perturbation  surface  fluxes  and  vertical  mixing  also  makes  the 
adjoint  quit  simple.  The  TLM  dependent  variables  appear  only  in  the  tri-diagonal  linear 
system  (82).  The  adjoint  is  therefore 

Ck-iXk-i  +  b^Xh  +  Uk-if-iXk+i  =  Xk,  (106) 

where  X^  is  now  the  vertical  profile  of  the  adjoint  variable  X  before  the  adjoint  mixing  is 
applied.  The  coefficients  are  identical  to  those  of  (83)  -  (89). 
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4.3  Vorticity  and  divergence  adjoint  equations 


We  need  the  adjoints  of  (75)  and  (76).  In  the  TLM  the  sequence  of  operations  is  to  evaluate 
(72),  (73),  and  (74)  in  grid  point  space,  add  the  diabatic  forcing  terms  and  perturbation 
geopotential  (j)j  and  then  do  grid-to-spectral  transforms  to  get  spectral  tendencies  of  (  and 
D.  In  the  adjoint,  therefore,  we  must  first  do  the  adjoint  of  the  transform.  As  mentioned 

earlier,  we  shall  postpone  the  details  of  this  for  now,  except  to  observe  that  a  grid-to- 

spectral  transform  adjoint  is  a  spectral-to-grid  transform  with  scaling.  The  reverse  is  true 
for  a  spectral-to-grid  transform  adjoint.  For  convenience  we  define  the  adjoint  operators 
Tta  and  «Sa,  where 

R  =  n^S  (107) 

for  the  grid-to-spectral  transform  adjoint,  and 

S  =  S^R  (108) 


for  the  spectral-to-grid  transform  adjoint. 

The  first  terms  to  consider  are  the  Laplacian  of  energy  terms  in  (74).  In  spectral 
space  we  have  an  adjoint  expression 

Sk  =  Sk  +  {dDk/dt) ,  (109) 


where  dDk/dt  corresponds  to  X'  in  (105).  Then  we  transform  Sk  to  grid  point  space  with 
the  appropriate  adjoint  transform  operator 


and  finally  the  adjoint  is 


Rk  =  n^Sk,  (110) 


Uk 

-  r/  +  R 

1  — 

(111) 

Vk 

1  —  11^ 

(112) 

<t>k 

=  (/>fc  +  Rki 

(113) 

Rk 

=  0, 

(114) 
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for  1  <  A:  <  L. 

Now  we  must  find  the  adjoints  of  a  HI.)  and  a  {HI.,  G'j.)  in  (75)  and  (76).  The  oc 
operator  is  a  spectral  transform  of  grid  point  vector  wind  components  to  spectral  vorticity 
and  divergence.  It  is  a  generalization  of  the  scalar  transforms  already  discussed  above 
and  will  be  described  in  greater  detail  in  a  section  below.  For  the  sake  of  brevity  we  will 
proceed  with  the  assumption  that  the  adjoint  of  the  a  operator  has  transformed  spectral 
adjoint  variables  {dC,k/dt)  and  {dDk/dt)  to  grid  point  Gk  and  Hk-  Then  the  adjoints  of  (72) 
and  (73)  are,  iox  \  <k  <  L\ 
vorticity 

Cfc  =  Ck  UkGk  +  VkHk,  (llh) 


U  wind 


Uk  =  Uk  + 


'f]^]  - 

L'5'?1a:+1/2  L' 


Ifc+l/Z  L'g^Jfc-1/2 

2Apf^ 


Hk  +  iCk  +  f)Gk, 


rr  tt  ,  L'9’7jfc-l/2  „■ 

c4-i  —  Gk-i  H  ttt::  71^, 

2Apfc 


Uk+i  —  74+1  “ 


A:-fl 


^  j  fc+1/2  rj- 


‘^^Pk 


V  wind 


Vk  =  Vk- 


ri^ 

l'^h+i/2 


~  b^]k- 


1/2 


2Apk 


Gk  +  (Cfc  +  f)Hk 


Vk-i  —  hfc_i  — 

14+1  =  14+1  + 


k-1/2 


“^^Pk 

W] 


fc+1/2 


Gk, 


Gk, 


vertical  velocity 


.  dp 

.  dp 

k-l/2 

,  Vk-Vk-i^  Uk-Uk-i„ 

T  „A_  brfc  A  —  -t/fc) 


k-1/2 


2Ap( 


2Api 


(116) 

(117) 

(118) 

(119) 

(120) 
(121) 

(122) 
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.  dp 

% 


J  fc+l/2 


.  dp 

% 


+ 


Vk+i-Vk^  Uk+i-Uk 


ik+l/2 


‘^^Pk 


-Gk  ~ 


2Apjfc 


Hk, 


and  if 


then  the  potential  temperature 


dk  —  + 


Exner  functions 


Cp 


Bk-l/2{Pk  —  Pk-112)  + -Bfc+i/2(Pfc+l/2  —  P  k) 

^Pk 


X\ 


(123)^ 


(124) 


Pk 

Pk-1/2 

Pk+1/2 


fc-1/2 


o?Xpk 

('pPk—lp^k 
a^Apk 


r.  ,  CpPk+lph^a 

+  a-^Ap,  ’ 


(125) 

(126) 
(127) 


terrain  pressure 


ABk 


TT  =  TT  — 


+ 


(Apfc)" 

.  dp 

k-1, 

.  dp 

-^^k  {Bk-lpiPk  —  Pk-112)  +  Bk+lpiP k+ll2  —  P k)) 


{Vk  -  Vk-,)  + 

I 

{Uk-Uk-i)  + 

'J  fc-1/2 

terrain  pressure  horizontal  gradients 


.  dp 


{Vk+i  -  Vk)  0.5Gfc 


fc+1/2 


.  dp 

^d-p 


{Uk+i  -  Uk)  0.5Hk 


k-\-l/2 


(128) 


diT  _  dn  Cp9k{l  —  p^)  Bk-ipjPk  ~  Pk-1/2)  +  Bk+i/2{P fc+1/2  —  Pk) 
dil  ~  djl^  ^2  [  Ap;- 


dn  _  dir  Cp9k 
dX  dX  a? 


Bk-i/2{Pk  —  P  fc-1/2)  +  ■^fc+1/2  (-Pfc+1/2  —  P  k) 

^Pk 


Hk. 


Gk,  (129) 
(130) 


Notice  how  the  adjoint  variables  are  accumulated  and  also  coupled  in  the  vertical. 
They  must  be  initialized  to  zero  at  the  beginning  of  each  adjoint  time  step  to  ensure  correct 
results.  In  the  actual  code  the  Uk,  Vk,  9k,  and  Qk  adjoints  are  initialized  from  the  solution 
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to  (106).  The  above  equations  are  the  result  of  the  adjoint  of  the  vorticity  and  divergence 
equations  only,  additional  sensitivities  will  be  accumulated  in  the  adjoint  equations  for  the 
other  dependent  variables.  The  actual  generation  of  the  adjoint  code  corresponding  to  these 
equations  was  by  the  Giering  and  Kaminski  (1996)  automatic  adjoint  software. 


4.4  Potential  temperature  adjoint  equations 

In  (70)  we  again  have  the  a  linear  operator  relating  spectral  temperature  tendencies  to 
horizontal  advection.  As  in  the  previous  section  we  assume  we  have  the  adjoint  of  this 
linear  operator,  so  the  spectral  adjoint  quantity  dOk/dt  yields  the  zonal  advection  adjoint 
Qk  and  the  meridional  advection  adjoint  h^-  Then  we  have  for  1  <  A:  <  L 


Ok  =  Ok  +  UkQk  +  y  khk,  (131) 

Uk  =  Uk  +  OkQki  (132) 

14  =  Vk  +  Okhk.  (133) 

Next  as  in  (110)  we  transform  dOk/dt  to  grid  point  space. 

Xk^n^idOk/dt),  (134) 


and  the  adjoint  equations  are; 
potential  temperature 


Dk 


Ok  =  Ok  + 


V 


§2 


■dp 

7]-^ 


^Pk 


4+1/2 


4-1/2 


=  0^ 


0^] 


k+1/2 


=  Ok-l/2  + 


k+1/2 


^Pk 


■q§£. 


1.- 


^^Jfc-1/2 

^Pk 


Xk  1 


Xk, 


(135) 

(136) 

(137) 
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vertical  velocity 


divergence 


■^P  _  ,  ^k-l/2-Ok 

9Pik-l/2  I  9p\k-l/2  ^Pk 

’.dp]  _  [.^pI  ,  ^k-dk^in 

^P\k+ll2  L  ^^Jfc+1/2  ^Pk 


Xk, 

Xk, 


Dk  —  Dk,  +  Ok  Xk, 


(138) 

(139) 

(140) 


terrain  pressure 

7^  ^  tt  +  ,^5\2  I  P^  {Ok-i/2-Ok)+  (0k-0k+ii2)\  Xk.  (141) 

{Xpk)  VL  ^P\k-ii2  L  / 

We  are  accumulating  sensitivity  contributions  into  the  adjoint  variables  to  add  to 

those  from  the  vorticity  and  divergence  adjoint  equations  of  the  previous  section.  Since  the 

moisture  adjoint  equation  for  q  is  exactly  analogous  to  the  6  equation  we  will  not  show  this 

set  of  equations.  The  model  code  for  them  is  available  upon  request  from  the  author. 


4.5  Hydrostatic  equation  adjoint 

The  adjoint  equations  of  the  hydrostatic  equation  (68)  are 


4>k+l 

—  <l>k+l  +  4'k, 

(142) 

4>k 

II 

(143) 

Ok 

—  Ok  +  {Pk+l/2  ~  P k)  <l>k, 

(144) 

Pk+ll2 

=  Pk+l/2  +  Ok  <l>k, 

(145) 

Pk 

~  Pk  Ok  (jki 

(146) 

for  1  <  A:  <  L  —  1  and 

Ok  ~  Ok  Cp(^Pk+l/2  Pk)  0fc) 


(147) 


Pk+l/2  —  Pk+1/2  +  Cpdk  4>k, 

(148) 

Pk  ~  Pk  Cp6k  (pki 

(149) 

for  k  =  L. 

Now  the  adjoints  of  the  vertical  temperature  interpolation  equations  (66)  and  (67): 


0  —ft  I  Pk-l/2  ^ 

—  f^k+l  T - ^  ^  f^k+1/2, 


P 


^c+l 


Pk 


f\  _  ft  \  ^k  n 

“k  —  “kP  ^k+\l2, 

Pk+l  —  Pk 

Pk  =  w> 

Pk+\  =  -ffc+1  ~  {^k  -  ^fc+l)  -p^Y  ^fc+1/2) 

^fc+1/2  =  0> 


(150) 

(151) 

(152) 

(153) 

(154) 


for  1  <  /c  <  L  —  1. 


4.6  Vertical  velocity  and  terrain  pressure  adjoints 

We  begin  the  calculations  for  the  adjoints  of  (61)  and  (62)  as  in  (110)  with  the  transform 
of  dr  jdt  to  grid  point  space,  combined  with  summations  over  the  previously  accumulated 
vertical  velocity  adjoint  variable.  We  have 


X, 


fc+l/2 


=  E 


dy 


n 


1  L  '^'/Jfc+l/2 


,  l<k<L-l, 


L-\ 


Xl+iI2  =  {dr/dt)  -  Bjt+l/2 

fc=l 


.  dp 

% 


(155) 

(156) 


J  fc+l/2 


Then  the  adjoint  expressions  are: 


Uk 


Vk 


ABk  dw 
. .  ^  dw 

Vk  —  AJSfc— Xfc_(-i/2, 


(157) 

(158) 
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Dk  =  Dk  —  Api^Xk+i/2, 

TT  =  TT  —  ABkDkXk+l/2-, 
dTT  _  dTt  ABk  YJ 

dX  ~  dx 

I  = 

-^fc-1/2  =  ■^fc-1/2  +  ^*+1/2) 


(159) 

(160) 
(161) 

(162) 

(163) 


ior  L  >  k  >  1.  Because  of  the  recursive  behavior  of  (163)  it  is  essential  to  evaluate  (158)- 
(163)  from  the  bottom  up,  i.e.  indexing  k  from  L  to  1. 


4.7  Exner  function  adjoints 

The  Exner  functions  from  (63)  and  (64)  are  dependent  only  on  the  terrain  pressure,  so 
the  adjoint  expressions  are  surprisingly  simple.  We  continue  accumulating  terrain  pressure 
sensitivity  with: 


,  Bk+i/2{Pk+l/2  —  Pk)  +  Bk-1/2{P k  —  P k-1/2)  c) 

TT  =  TT  -I - - - T-r - -nb 

_  ^Pk 

+  (164) 

Pk+ll2 

for  1  <  k  <  L. 

At  this  point  we  have  completed  grid  point  space  computations  of  all  adjoint  depend¬ 
ent  variables.  In  summary,  since  the  TLM  completes  a  time  step  with  the  leapfrog  and  time 
filter  operations,  we  began  the  process  with  the  adjoint  of  these  time  stepping  algorithms 
with  adjoint  variables  in  spectral  space.  We  then  transformed  these  to  grid  point  space 
using  a  spectral-to-grid  point  adjoint  transform,  followed  be  the  adjoint  calculations  de¬ 
scribed  above  in  exact  reverse  sequence  from  the  TLM  order  of  calculations.  The  next  step 
is  to  transform  back  to  spectral  space  using  a  grid-to-spectral  adjoint  transform. 
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4.8  Spectral  transform  adjoints 

The  scalar  grid  point  to  spectral  transform  (41)  can  be  written  symbolically  in  matrix  form 


as 


5-1 

Sf/'i 

(165) 


where  K  is  the  linear  operator  composed  of  the  discrete  Fourier  transforms  and  the  Gaussian 
quadratures  which  converts  a  horizontal  grid  point  field  of  mn  degrees  of  freedom  to  a  set 
of  complex  spherical  harmonic  coefficients  with  degrees  of  freedom.  From  the  definition 
of  the  adjoint  we  have 


All 

^mn 

(166) 


where  the  (  )  are  now  adjoint  variables  and  K'^  is  the  matrix  transpose  of  K.  After  some 
tedious  but  straightforward  algebra  (166)  can  be  shown  to  be  equivalent  to 

^  M 

m=—M 

Equation  (166)  is  identical  in  form  to  (40)  except  for  the  Gaussian  weights  Wj  proceeding  the 
summations.  Therefore  the  adjoint  of  the  grid-to-spectral  transform  is  the  spectral-to-grid 
transform  of  the  forward  model  with  an  additional  scaling  factor  weighting  each  Gaussian 
latitude  of  the  horizontal  grid  point  field  after  the  transform. 

A  similar  examination  yields  the  adjoint  of  the  spectral-to-grid  transform  as: 

(3M+l)/2  _ 

s^=  Y.  (168) 

which  is  similar  to  (43)  except  that  the  Gaussian  quadrature  weights  are  absent  from  the 
summation.  In  the  adjoint  code  the  grid-to-spectral  transform  code  of  the  forward  model  is 


M 

E 

|n=|m| 


5771  ryn 

n 


\n) 


AmXi 


(167) 


36 


used  except  that  an  array  with  element  values  of  1.0  instead  of  the  quadrature  weights  is 
passed  as  a  subroutine  argument. 

Not  surprisingly,  the  adjoint  of  the  spectral  vorticity  and  divergence  to  grid  point 
wind  transform  is  a  grid-to-spectral  transform:  analogous  to  (52)  and  (53).  Deriving  the 
adjoint  expressions  in  matrix  form  as  shown  above,  it  can  be  shown  that 

1 

n{n  +  1) 


c 


m 

n 


(3M-fl)/2 

Y.  {im^"  [r{M,)]  ?„"■(«) 


+  (i- 


and 


=  - - - 

n(n  +  1) 


(3M+l)/2 

X:  [t/(M,)] 

j~l 


(169) 


(170) 


Similarly,  the  adjoint  of  the  grid  point  wind  to  spectral  vorticity  and  divergence  is  a 
spectral-to-grid  transform  analogous  to  (50)  and  (51)  and  is 


U {Xii  fij')  —  Wj 
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M 
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m^—M  |_n=:|m| 
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'j  m^  —  M 
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(171) 


and 


V{Xi,iJ,j)  =  Wj 


1-fj, 


M 
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-im 


'3  m=—M 


M 


Y  CPTM 

n=\m\ 


^imXi 


M 

-  E 

m=—M 


^  ^  dP'^(a  ) 

n=imj 


dji 


J.m\i 


(172) 


Comparing  (169)  and  (170)  to  (52),  (53),  and  (43)  suggests  we  scale  the  input  U  and 
V  with  1  —  }Xj,  set  the  quadrature  weights  to  1.0,  call  the  forward  model  subroutine  for 
grid  point  wind  to  spectral  vorticity  and  divergence  transforms,  and  then  scale  the  output 
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spectral  coefficients  with  1.0/(n(n  +  1))  to  produce  the  correct  adjoint  values  of  C  and  D. 
Likewise  for  (171)  and  (171)  we  scale  the  input  C  and  D  with  n{n  +  1),  pass  these  to  the 
forward  model  spectral  vorticity  and  divergence  to  grid  point  wind  transform  subroutine, 
and  then  scale  the  output  grid  point  winds  with  Wj/{1  —  fXj)  to  yield  the  correct  adjoint 
values  U  and  V.  Note  the  symmetry  of  the  scaling  factors  between  the  two  transform 
directions.  This  is  necessary  if  the  adjoint  transforms  are  to  retain  the  inverse  properties 
of  the  companion  transforms  of  the  forward  model. 


5  Testing  of  TLM  and  adjoint  models 

Developing  tangent-linear  and  adjoint  models  from  a  parent  non-linear  NWP  model 
such  as  NOGAPS  requires  an  orderly  sequence  of  steps  to  minimize  the  risk  of  errors.  For¬ 
mulation  of  the  TLM  model  described  in  section  3  is  straightforward  but  tedious.  Recently 
automatic  TLM  generation  software  has  been  introduced  (Giering  and  Kaminski,  1996), 
but  it  was  not  available  for  the  NOGAPS  TLM  ,  so  all  coding  was  by  hand.  The  testing 
method  for  NOGAPS  TLM,  as  with  any  linearized  model,  is  to  examine  its  behavior  for  a 
range  of  perturbation  initial  conditions.  The  criteria  used  is  that 

Le  ^  N{x  +  e)  —  N{x)  (173) 

where  L  is  the  tangent  linear  operator  derived  from  the  non-linear  model  N  and  e  is  the  initial 
perturbation  on  a  full  initial  field  x.  In  the  limit  of  sufficiently  small  e,  (173)  will  become 
an  effective  equality  if  the  linearization  is  robust.  As  e  increases  the  non-linearity  of  N  will 
cause  the  approximation  to  break  down.  Errico  et  al.(1993)  shows  that  for  an  e  equivalent 
to  analysis  uncertainty  (errors),  the  TLM  of  a  mesoscale  forecast  model  (MM4)  is  quite 
accurate  for  48  and  even  72  hours  for  winter  forecast  cases.  We  expect  the  NOGAPS  TLM 
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TLM  ERROR 


Figure  2:  Correlation  of  TLM  solution  to  difference  between  two  non-linear  solutions  for  12 
hour  forecast  lengths.  Model  resolution  is  T21. 

to  have  similar  capabilities.  To  test  the  NOGAPS  TLM  we  perturb  each  spectral  coefficient 
of  each  dependent  variable  randomly  with  values  in  the  range  ±e.  The  correlation  error 
between  the  NOGAPS  TLM  and  difference  of  two  non-linear  (dry)  forecasts  as  a  function 
of  e  is  shown  in  Fig  2. 

Fig  2  is  a  log-log  plot  of  (1 -correlation  coefficient)  vs.  perturbation  size.  Over  most 
of  the  perturbation  size  range  the  relationship  is  quite  linear,  corresponding  to  a  quadratic 
growth  of  the  correlation  error  as  a  function  of  e.  In  fact,  it  shows  that  the  non-linear 
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solutions  of  the  dry  dynamics  have  essentially  quadratic  behavior  for  perturbations  with 
magnitudes  nearly  as  large  as  the  original  unperturbed  coefficients.  It  also  shows  that 
the  TLM  solution  is  quite  accurate  over  this  range,  with  better  than  95%  correlation  over 
the  same  perturbation  range.  Only  for  extremely  large  perturbations  does  the  quadratic 
behavior  break  down,  probably  because  the  cubic  dependence  of  the  vertical  advection  terms 
begins  to  have  significant  impact  on  the  non-linear  solutions.  Certainly  the  lack  of  any  moist 
physics  or  other  discrete  physical  processes  contributes  to  this  surprising  result,  and  if  we 
repeated  the  experiment  for  a  difference  between  two  full  physics  non-linear  forecasts  we 
would  certainly  see  a  breakdown  of  the  TLM  solutions  much  sooner. 

Adjoint  model  testing  presents  special  problems  because,  unlike  forward  model  testing, 
there  is  no  simplified  physical  or  meteorological  solutions  available  that  allow  validation  of 
some  basic  model  properties  such  as  dynamic  balance  or  energy  conservation.  Therefore  we 
must  depend  on  linear  algebra  tests  derived  from  the  formal  transpose  relationship  between 
the  TLM  model  and  its  companion  adjoint  model.  Insuring  that  the  TLM  is  a  correct  and 
robust  proxy  for  the  parent  nonlinear  model,  using  the  method  described  above,  is  also  a 
top  priority  for  developing  a  useful  adjoint  model. 

The  diagram  below  shows  schematically  a  forward  (TLM)  integration  from  time  0 
to  time  t,  and  a  backward  (ADJ)  adjoint  integration  from  time  t  to  time  0.  The  adjoint 
integration  is  shown  in  terms  of  a  gradient  dJ/dX  of  a  cost  function  J  whose  gradients  are 
of  interest.  Typically  J  is  based  on  an  energy  norm  or  some  integrated  measure  of  forecast 
skill  such  as  RMS  error.  See  Ehrendorfer  and  Errico  (1995)  for  a  discussion  of  the  choice 
of  norms. 


Xo 


TLM 


dJ 

dXO 


ADJ 


dJ 

axt 


In  this  diagram  the  AT  and  dJ/dX  represent  vectors  of  spectral  coefficients.  For  an 
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adjoint  model  to  be  the  correct  transpose  of  a  TLM  model,  we  must  have 

where  <  >  represents  the  inner  product.  This  is  called  the  gradient  test,  and  is 

performed  by  making  the  TLM  and  adjoint  integrations  shown,  saving  the  initial  and  final 
conditions  as  indicated  so  the  inner  products  can  be  computed.  For  testing  purposes  there 
is  no  constraint  on  the  nature  of  the  Xq  and  dJ/dXt  initial  conditions,  for  the  NOGAPS 
adjoint  model  random  numbers  were  used. 

The  gradient  test  is  an  extremely  sensitive  indicator  of  the  validity  of  an  adjoint  code. 
For  64  bit  floating  point  calculations  the  two  inner  products  above  should  have  9-10  decimal 
digit  agreement  for  real  data  integrations  of  24-48  hours.  Poorer  agreement  than  this  is 
almost  certainly  due  to  coding  errors,  often  because  accumulated  adjoint  sensitivities  are 
not  initialized  properly.  In  practice  it  is  usually  best  to  first  do  gradient  testing  of  individual 
code  modules  for  single  time  steps,  combining  modules  as  they  are  validated.  When  the 
complete  adjoint  model  is  assembled  and  tested  for  a  single  time  step,  then  multiple  time 
step  testing  can  begin.  This  is  the  best  way  to  isolate  errors  in  variable  initialization. 

6  The  singular  vector  system 

A  powerful  application  of  TLM  and  adjoint  models  is  their  combination  to  form  a 
singular  vector  matrix  operator,  the  leading  eigenvalues  and  eigenvectors  of  which  describe 
the  most  unstable  dynamical  structures  of  a  time-evolving  basic  state.  This  instability  is 
manifested  by  the  growth  of  some  appropriately  chosen  norm  computed  from  the  perturb¬ 
ation  variables  of  the  TLM  and  adjoint  models.  The  total  energy  norm  is  most  commonly 
used.  It  is  the  sum  of  perturbation  kinetic  energy  and  available  potential  energy  summed 
over  all  resolved  wavenumbers  and  the  depth  of  the  atmosphere. 
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6.1  Energy  norm  formulation 


At  any  time  t  the  energy  norm  is 

M  M  r  L 


where 


m=-M  n=m  ^k=l 

+  a{k)e^{k,ty  +/?7r;^(t)2| 


(175) 
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with  5  =  4  for  m  =  0  and  5  =  2  for  m  >  0.  Po{k),To,  and  ttq  are  suitable  chosen  reference 
values  for  Exner  function,  temperature,  and  terrain  pressure,  respectively. 

In  compact  matrix  form  (175)  is 


Et  =  xfC^Xt, 


(176) 


where  the  x  is  a  vector  of  all  spectral  coefficient  values  of  the  perturbation  variables  C:D,9, 
and  7r,  and  is  a  diagonal  matrix  of  the  energy  norm  coefficients  7,  a,  and  p.  At  time 
t  =  0  the  energy  is 


Eo  =  xfC'^xo. 

From  the  definitions  of  the  TLM  and  adjoint  model  we  have 

(177) 

Xt  =  Lxo, 

(178) 

xj  =  xfL^, 

(179) 

which  can  be  substituted  into  (176)  to  yield  the  energy  att  =  t. 

Et  =  x^L^C^Lxo. 


(180) 
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6.2  Singular  vector  formulation 

We  want  to  find  the  eigenvalues  and  eigenvectors  of  the  linear  system  (180)  which  maximize 
the  growth  of  E  over  the  time  interval  t.  If  we  normalize  the  initial  vector  xq  such  that 


Eo  =  1,  (181) 

then  we  can  pose  the  algebraic  eigenvalue  problem  in  variational  form  as 

Jt  =  x^L^C^Lxo  -  Xix^C^xo  -  1).  (182) 

(183) 
However, 

Xo  =  C"^2:o,  (184) 

and  multiply  both  sides  by  we  get 

C'^L^C^LC-^zo  -  Azo  =  0,  (185) 

a  much  more  tractable  symmetric  eigenvalue  problem.  Notice  that  substitution  of  (184)  into 
(177)  and  imposing  (181)  yields 

2:^20  =  1-  (186) 


Taking  the  first  variation  of  (182)  yields 

=  L^C^Lxo  -  AC^xo  =  0. 

OXo 

In  this  form  (183)  is  awkward  because  it  is  not  a  symmetric  eigenvalue  problem, 
if  we  introduce  an  energy  norm  scaling  transformation 


6.3  The  Lanczos  algorithm 

The  primary  obstacle  to  solution  of  (185)  is  the  enormous  size  of  the  matrix  operator.  As 
described  in  section  1,  it  can  easily  be  of  0(10’’)  for  linear  operators  based  on  models  such 
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as  NOGAPS.  However,  because  we  are  interested  only  in  a  few  of  the  leading  eigenvalues 
of  the  system,  we  do  not  need  an  explicit  representation  of  the  operator.  The  Lanczos 
algorithm  (see  Golub  and  Van  Loan  (1989),  Strang  (1986),  and  Simon  (1984)  for  details)  is 
ideally  suited  for  this  problem,  since  it  only  requires  as  its  input  the  result  vector  produced 
when  a  normalized  input  vector  zq  is  passed  through  the  matrix  operator  in  (185).  For 
computation  purposes  the  matrix  operator  is  treated  as  a  ’black  box  called  by  the  Lanczos 
algorithm  code  as  it  iterates. 

The  Lanczos  iteration  process  proceeds  as  follows,  given  an  initial  zq  scaled  to  satisfy 

(186); 

1.  Compute  C“^zo-  This  is  equivalent  to  scaling  the  zq  back  to  the  model  dependent 
variables  (,  D,  9,  and  tt. 

2.  Compute  LC^^^o-  This  is  the  TLM  integration  over  an  appropriate  forecast  time, 
called  the  optimization  time. 

3.  Compute  C^LC'^^zo.  This  is  the  energy  norm  scaling. 

4.  Compute  L^C^LC'^zq.  This  is  the  adjoint  model  integration  over  the 
optimization  time  . 

5.  Compute  C'^L^C^LC'^zo.  This  transforms  the  dependent  variables  back 
to  the  zo  energy  norm  scaling  form. 

The  result  of  step  5  is  passed  to  the  Lanczos  algorithm,  which  re-scales  the  vector  to 
yield  an  updated  estimate  for  zo,  and  then  steps  1  -  5  are  repeated.  The  Lanczos  algorithm 
generates  estimates  for  the  largest  eigenvalues  and  associated  eigenvectors  as  these  iterations 
proceed,  both  refining  the  estimates  and  increasing  the  number  of  candidate  eigenvalues. 
Typically  the  number  of  converged  eigenvalues  equals  about  1/3  the  number  of  iterations. 
The  Lanczos  calculations  themselves  are  not  computationally  expensive,  but  since  each  it¬ 
eration  involves  both  a  TLM  and  adjoint  model  time  integration,  the  overall  computational 
cost  is  considerable.  For  example,  A  50  iteration  singular  vector  T47  NOGAPS  calculation 
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for  a  36  hour  optimization  time  requires  about  25000  CPU  seconds  on  a  CRAY  COO  super¬ 
computer.  This  includes  the  cost  of  a  non-linear  forecast  to  produce  the  time  dependent 
basic  state  trajectory.  I/O  costs  are  also  large,  since  the  TLM  and  adjoint  model  each  must 
read  a  ~  100  Mbyte  trajectory  file  every  iteration. 

6.4  The  local  projection  operator 

As  described,  the  singular  vector  system  implicitly  defines  a  global  optimization  domain 
because  of  the  global  nature  of  the  spherical  harmonic  coefficients.  More  desirable  is  a 
regionally  defined  optimization  domain  to  allow  exclusion  of  singular  vectors  which  lie  out¬ 
side  an  area  of  interest.  Such  a  targeted  domain  will  be  influenced  by  a  small  number  of 
singular  vectors,  typically  1-3,  requiring  only  about  10  Lanczos  iterations,  instead  of  the 
~  50  iterations  needed  to  find  the  leading  singular  vectors  for  a  global  domain. 

A  regional  optimization  domain  is  constructed  by  simply  discarding  the  variance  in 
the  model’s  dependent  variables  outside  the  domain,  so  that  energy  norm  growth  will  depend 
only  on  what  happens  inside  the  domain.  We  define  a  local  projection  operator  (LPO)  as 

F  =  S-^nS,  (187) 

where  S  is  the  spectral-to-grid  transform  operator  (40),  its  inverse  (41),  and  77.  is  a 
Gaussian  grid  point  mask  equal  to  1.0  in  the  optimization  domain  and  0.0  elsewhere.  The 
mask  77  is  a  grid  point  operator,  but  the  LPO  is  an  operator  in  spherical  harmonic  space, 
so  in  practice  there  is  some  non-zero  grid  point  variance  due  to  spectral  truncation  Gibb’s 
phenomena  in  the  periphery  of  the  projection  area.  This  is  only  cosmetic,  since  the  variance 
drops  off  very  rapidly  away  from  the  projection  area  and  does  not  contribute  significantly 
to  perturbation  energy  growth.  F  is  a  symmetric  operator,  so  it  can  be  introduced  into 
(185),  yielding  the  modified  symmetric  eigenvalue  problem 

C'^L^FC^FLC-^zo  -  Azo  =  0,  (188) 
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that  is  still  solvable  with  the  Lanczos  algorithm. 


7  Example  Calculations 

We  will  demonstrate  the  application  of  the  TLM,  adjoint,  and  singular  vector  model¬ 
ing  systems  on  a  case  study  from  the  recent  Fronts  &  Atlantic  Storm  Tracks  Experiment 
(FASTEX)  (Joly  et  al  1997  ).  The  particular  example  includes  the  ’FASTEX’  cyclone,  an 
intense  and  rapidly  moving  storm  that  developed  in  the  mid-Atlantic  between  1997021612 
UTC  and  1997021912  UTC  and  was  the  subject  of  intense  observational  study  during  its 
life  cycle.  Figures  3a,  3b,  3c,  3d,  3e,  and  3f  are  the  NOGAPS  1997021612  UTC  initial  con¬ 
ditions  and  the  24,  36,  48,  60,  and  72  hour  forecasts  respectively.  The  incipient  FASTEX 
cyclone  is  initially  visible  at  tau=24  as  a  diffuse  low  pressure  area  south  of  the  Canadian 
Maritime  provinces.  During  the  next  48  hours  the  storm  races  across  the  Atlantic  and  be¬ 
comes  a  relatively  small,  very  intense  system  that  brought  gale  force  winds  to  Ireland,  Great 
Britain,  and  the  North  Sea.  This  storm  is  ideally  suited  for  study  with  singular  vectors  and 
adjoint  sensitivities  because  of  its  relatively  small  size  in  the  FASTEX  validation  area  of 
the  eastern  Atlantic.  The  72  hour  1997021612  UTC  forecast  shown  here  is  the  basic  state 
trajectory  chosen  for  the  singular  vector  calculations. 

7.1  Singular  vectors 

During  FASTEX  singular  vector  and  adjoint  sensitivity  calculations  were  used  to  select 
areas  for  objective  targeting  of  aircraft  dropsondes  in  the  far-upstream  area  from  the  val¬ 
idation  zone.  Because  of  aircraft  scheduling  constraints  and  experimental  planning  require¬ 
ments  it  was  necessary  to  run  the  singular  vector  calculations  more  than  24  hours  before 
the  actual  targeting  times.  In  the  example  shown  here  a  36  hour  lead  time  was  used,  giving 
a  36  hour  optimization  time  for  the  singular  vectors  from  1997021800  UTC  to  1997021912 
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Figure  3:  Time  series  of  sea  level  pressure  forecasts  showing  72  hour  development  history 
of  the  ’FASTEX’  cyclone. 


UTC  of  the  trajectory  forecast.  The  ’target  time’  for  the  aircraft  mission  was  therefore 
1997021800  UTC  (Tau=36  in  the  trajectory  forecast),  and  the  verification  time  1997021912 
UTC  (Tau=72).  Figure  4a  is  the  perturbation  temperature  variable  of  the  leading  singular 
vector  at  680  mbs  for  the  LPO  box  shown  bounding  the  FASTEX  validation  area.  It  has 
the  largest  growth  of  the  total  perturbation  energy  norm  in  the  LPO.  Figure  4b  is  this  per¬ 
turbation  temperature  after  a  36  hour  integration  of  the  TLM  initialized  with  the  leading 
singular  vector.  It  is  the  ’evolved’  singular  vector,  and  shows  the  growth  of  the  perturbation 
energy  in  the  LPO. 

SINGULAR  VECTOR  at  TARGETED  OBSERVING  TIME  SINGULAR  VECTOR  at  FORECAST  VERIFICATION  TIME 


Figure  4:  (a)  Leading  singular  vector,  and  (b)  its  36  hour  TLM  evolution  for  the  FASTEX 
LPO.  Note  that  the  contour  interval  for  the  evolved  vector  is  10  times  that  for  the  initial 
vector. 

The  potential  of  singular  vector  calculations  such  as  this  is  considerable.  They  tell 
us  ’where  the  action  is’  in  our  forecast  models,  and  by  implication,  in  the  real  atmosphere. 
The  use  of  an  LPO  allows  us  to  focus  our  attention  on  a  forecast  area  of  interest.  NWP 
model  forecast  error  is  the  result  of  initial  condition  error  growing  during  a  forecast.  This 
error  growth  is  concentrated  in  areas  of  atmospheric  instability  such  as  baroclinic  zones. 
The  singular  vectors  tell  us  where  these  areas  of  instability  are,  and  evolved  vectors  can 
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tell  us  where  the  forecast  error  will  occur.  Because  the  singular  vectors  are  normalized 
eigenvectors,  we  cannot  infer  any  magnitude  or  sign  information  about  the  forecast  error, 
at  least  in  real  time,  but  they  do  tell  us  where  improving  our  initial  conditions,  i.e.,  reducing 
the  initial  analysis  error,  will  have  the  greatest  potential  for  reducing  forecast  error.  This 
was  the  premise  of  the  objective  targeting  strategy  in  the  far-upstream  area  of  FASTEX, 
where  aircraft  dropsondes  were  targeted  in  the  areas  of  maximum  singular  vector  amplitude. 


7.2  Adjoint  sensitivities 


Another  equally  powerful  application  of  the  adjoint  model  is  computing  the  sensitivities  of 
some  appropriately  chosen  cost  function  J  to  perturbations  in  meteorological  variables.  In 
this  example  we  define  a  cost  function  as  the  perturbation  vorticity  in  the  lower  troposhere 
of  the  FASTEX  LPO  at  a  time  of  interest,  e.g.  the  verification  time  (1997021912  UTC)  of 
our  example  NWP  forecast.  Figure  5a  is  this  cost  function.  We  ask  the  question:  what  will 
be  the  sensitivity  dJ/dTo  of  this  cost  function  to  perturbations  to  an  atmospheric  variable, 
in  this  case  temperature,  at  some  earlier  time,  i.e.,  at  1997021800  UTC?  The  adjoint  model 
integration  is 


dJ  _  y 

Wt’ 


(189) 


where  the  ’initial’  condition  dJ/dTt,  is  the  derivative  of  our  perturbation  vorticity  cost 
function  at  the  verification  time  1997021912.  The  sensitivity  pattern  shown  in  Figure  5b 
is  dJ/dTo  at  680  mbs.  Because  we  use  the  same  basic  state  trajectory  as  in  the  singular 
vector  calculations,  and  the  cost  function  is  non-zero  over  the  same  LPO  used  for  the 
singular  vector  calculations,  the  resulting  patterns  are  strikingly  similar  to  the  singular 
vector  temperatures,  although  of  larger  horizontal  scale.  Adjoint  sensitivity  patterns  such 
as  this  can  be  interpreted  as  a  composite  of  the  full  spectrum  of  singular  vectors,  essentially 
an  weighted  average  depending  on  the  eigenvalues  of  each  singular  vector,  which  tends  to 
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spread  them  out  as  we  see.  The  choice  of  cost  function  is  obviously  an  important  factor  in 
this  interpretation,  but  that  is  beyond  the  scope  of  this  report. 


INITIAL  ADJOINT  FIELD  at  FORECAST  VERIFYING  TIME  ADJOINT  SENSITIVITY  AT  TARGETED  OBSERVING  TIME 


COST  FUNC:  Initial  VOR  680hPa  d(J)/d(TEMP)  680hPa 


Figure  5:  (a)  Perturbation  vorticity  cost  function  J  in  FASTEX  LPO  at  1997021912  UTC, 
(b)  dJ/dT  at  1997021800  UTC. 


8  Summary 

The  NOGAPS  adjoint  modeling  system  is  a  powerful  new  tool  for  understanding  the  beha¬ 
vior  of  the  atmosphere  and  the  NWP  models  used  to  predict  it.  In  recent  years  the  major 
operational  NWP  centers  have  reached  a  plateau  in  forecast  model  skill,  particularly  at 
forecast  times  less  than  72  hours.  Increasing  model  resolution  with  each  new  generation 
of  computer  technology,  once  a  near  guarantee  of  better  model  performance,  no  longer  has 
the  positive  impact  it  once  did.  There  are  many  reasons  for  this  plateau  in  NWP  model 
skill,  including  unrealistic  parameterization  of  diabatic  processes,  model  climate  drift  that 
is  insensitive  to  resolution  change,  and  inadequate  date  coverage  in  the  tropics  and  southern 
hemisphere.  Short-term  forecast  error  (<  72  hours),  however,  is  primarily  an  initial  value 
problem,  and  reducing  initial  analysis  error  is  the  only  way  to  improve  these  forecasts. 


Singular  vector  and  adjoint  sensitivity  calculations  such  as  the  examples  in  section  7  can 
now  find  the  areas  of  collocation  of  analysis  error  and  sensitivity  to  these  errors.  We  can 
focus  our  attention  on  these  ’hot  spots’  with  additional  observations,  special  attention  to 
data  quality  control,  or  special  objective  analysis  techniques,  to  reduce  the  analysis  error. 

An  exhaustive  discussion  of  other  potential  applications  of  adjoint  models  is  far  beyond 
the  scope  of  this  report.  Errico  (1997)  gives  very  complete  and  detailed  descriptions  of  the 
many  other  ways  adjoints  have  been  used  to  open  new  meteorological  research  areas  or 
re-address  old  questions  in  more  rigorous  and  efficient  ways.  This  report  should  give  a 
flavor  of  what  real  TLM  and  adjoint  models  are,  and  why  they  will  be  so  important  for  the 
future  of  meteorological  research. 
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